#Program for analysis of population growth rate
rm(list=ls())
library(boot);library(popbio);library(rjags);library(coda)
#root = "//Volumes/wcnr-network/Research/Hobbs/A_Bison_2011/Analysis for Ecological Monograph/"
root = "/Users/Tom/Documents/YNP_Bison_Model/"

path = function (stem,r=root){
return(paste(r,stem, sep=""))	
}

load(path("Model and analysis code/Frequency Dependent/One_beta_coda.Rdata"))
load(path("Model and analysis code/Frequency Dependent/One_beta_jags.Rdata"))

setwd(path("Model and analysis code/Analysis/R_0 and lambda/"))

set.seed(1)
#plotting code at bottom






#rename MCMC output and covert to dataframe

b=One_beta.coda
d=rbind(b[[1]],b[[2]],b[[3]])
c=as.data.frame(d)
n.rep=50000
#storage vectors and matrices
n2=n.rep
R0=numeric(n2)
p=numeric(2)
phi=numeric(1)
f.con=numeric(n2)
f.neg=numeric(n2)
f.pos=numeric(n2)
lambdaI=numeric(n2)
lambdaH=numeric(n2)
lambda.diff=numeric(n2)
adult.fem = numeric(n2)	
yearling.fem = numeric(n2)
adult.male = numeric(n2)
calf = numeric(n2)

adult.femI = numeric(n2)	
yearling.femI = numeric(n2)
adult.maleI = numeric(n2)
calfI = numeric(n2)


stage=matrix(nrow=n2,ncol=8)
stageI=matrix(nrow=n2,ncol=8)

b=numeric(3)
O=matrix(nrow=n2,ncol=15)



##########analysis of lambda for single sex populations
for(j in 1:n2){
	A=matrix(0,8,8)
	
	t=1

	#p=numeric(2)
	phi=numeric(1)
	b=numeric(4)
#make a random draw to index MCMC output
	i=round(runif(n=1, min=1, max= nrow(c)))

#random draw of parameter values
	beta=c$beta[i]
	#survivals
	p[1]=c$"p[1]"[i]
	p[2]=c$"p[2]"[i]
	p[3]=c$"p[3]"[i]
	psi=c$psi[i]
	v=c$v[i]
	#females only model
	m = 0
	f.pos=c$"b[2]"[i] *.5
	f.neg = c$"b[1]"[i]*.5
	f.con = c$"b[3]"[i]*.5
	
	
 #Average probability of infection
	phi[t]=One_beta.jags$xphi[i]
	
#Projection matrix
A[1,3]<-p[2]*f.neg[t]*(1-phi[t])
A[1,6]<-p[2]*f.con[t]*(1-v)*(1-phi[t])
A[1,7]<-p[2]*(f.pos[t] *(1-psi)*(1-phi[t])  +  f.pos[t] * psi*(1-v)*(1-phi[t]))


#row 2
A[2,1]<-p[1]*(1-m)*(1-phi[t])
#row 3
A[3,2]<-p[2]*(1-phi[t])
A[3,3]<-p[2]*(1-phi[t])

#row4 
A[4,3]<-p[2]*f.neg[t]*phi[t]
A[4,6]<-p[2]*f.con[t]*(v+phi[t]-v*phi[t])
A[4,7]<-p[2]*(f.pos[t]*psi*(v+phi[t]-v*phi[t]) + f.pos[t]*(1-psi)*phi[t])

#row5
A[5,1]<-p[1]*phi[t]*(1-m)
A[5,4]<-p[1]*(1-m)

#row6
A[6,2]<-p[2]*phi[t]
A[6,3]<-p[2]*phi[t]
A[6,5]<-p[2]
A[6,7] <- 0

#row 7
A[7,6]<-p[2]
A[7,7]<-p[2]


#row 8
A[8,1]<-p[1]*m
A[8,4]<-p[1]*m
A[8,8]<-p[3]

ea=eigen.analysis(A)

 	#Calculate growth rate for all female population
 	lambdaI[j]=ea$lambda1


#healthy populaiton
	f.pos=c$"b[1]"[i]*.5
	f.neg = c$"b[1]"[i]*.5
	f.con = c$"b[1]"[i]*.5
	m =0
	phi[t] = 0
	v = 0
	psi = 0	

A[1,3]<-p[2]*f.neg[t]*(1-phi[t])
A[1,6]<-p[2]*f.con[t]*(1-v)*(1-phi[t])
A[1,7]<-p[2]*(f.pos[t] *(1-psi)*(1-phi[t])  +  f.pos[t] * psi*(1-v)*(1-phi[t]))



#row 2
A[2,1]<-p[1]*(1-m)*(1-phi[t])
#row 3
A[3,2]<-p[2]*(1-phi[t])
A[3,3]<-p[2]*(1-phi[t])

#row4 
A[4,3]<-p[2]*f.neg[t]*phi[t]
A[4,6]<-p[2]*f.con[t]*(v+phi[t]-v*phi[t])
A[4,7]<-p[2]*(f.pos[t]*psi*(v+phi[t]-v*phi[t]) + f.pos[t]*(1-psi)*phi[t])


#row5
A[5,1]<-p[1]*phi[t]*(1-m)
A[5,4]<-p[1]*(1-m)

#row6
A[6,2]<-p[2]*phi[t]
A[6,3]<-p[2]*phi[t]
A[6,5]<-p[2]
A[6,7] <- 0

#row 7
A[7,6]<-p[2]
A[7,7]<-p[2]


#row 8
A[8,1]<-p[1]*m
A[8,4]<-p[1]*m
A[8,8]<-p[3]

	ea=eigen.analysis(A)
 	#Get growth rate for females only population
 	
 lambdaH[j]=ea$lambda1
lambda.diff[j] = lambdaH[j] - lambdaI[j]
	
}


##########analysis of stable stage distribution for two sex model
for(j in 1:n2){
	A=matrix(0,8,8)
	
	t=1

	#p=numeric(2)
	phi=numeric(1)
	b=numeric(4)

	i=round(runif(n=1, min=1, max= nrow(c)))

	
	beta=c$beta[i]
	
	p[1]=c$"p[1]"[i]
	p[2]=c$"p[2]"[i]
	p[3]=c$"p[3]"[i]
	psi=c$psi[i]
	v=c$v[i]
	
	m = c$m[i]
	f.pos=c$"b[2]"[i] 
	f.neg = c$"b[1]"[i]
	f.con = c$"b[3]"[i]
	
	
 
	phi[t]=One_beta.jags$xphi[i]
	

A[1,3]<-p[2]*f.neg[t]*(1-phi[t])
A[1,6]<-p[2]*f.con[t]*(1-v)*(1-phi[t])
A[1,7]<-p[2]*f.pos[t]*(1-psi)*(1-phi[t])  #changes


#row 2
A[2,1]<-p[1]*(1-m)*(1-phi[t])
#row 3
A[3,2]<-p[2]*(1-phi[t])
A[3,3]<-p[2]*(1-phi[t])

#row4 
A[4,3]<-p[2]*f.neg[t]*phi[t]
A[4,6]<-p[2]*f.con[t]*(v+phi[t]-v*phi[t])
A[4,7]<-p[2]*f.pos[t]*psi*(v+phi[t]-v*phi[t]) #

#row5
A[5,1]<-p[1]*phi[t]*(1-m)
A[5,4]<-p[1]*(1-m)

#row6
A[6,2]<-p[2]*phi[t]
A[6,3]<-p[2]*phi[t]
A[6,5]<-p[2]
A[6,7] <- 0

#row 7
A[7,6]<-p[2]
A[7,7]<-p[2]


#row 8
A[8,1]<-p[1]*m
A[8,4]<-p[1]*m
A[8,8]<-p[3]

ea=eigen.analysis(A)
 	#browser()
 	#Calculate growth rate for all female population
 	lambdaI[j]=ea$lambda1
	#browser()
	
	stage[j,]=ea$stable.stage
	
 	stageI[j,]=ea$stable.stage
	adult.femI[j] = sum(stageI[j,c(3,6,7)])
	yearling.femI[j] = sum(stageI[j,c(2,5)])
	adult.maleI[j] = stageI[j,8]
	calfI[j] = sum(stageI[j,c(1,4)])


#healthy populaiton
	f.pos=c$"b[1]"[i]
	f.neg = c$"b[1]"[i]
	f.con = c$"b[1]"[i]
	m = c$m[i]
	phi[t] = 0
	v = 0
	psi = 0	

A[1,3]<-p[2]*f.neg[t]*(1-phi[t])
A[1,6]<-p[2]*f.con[t]*(1-v)*(1-phi[t])
A[1,7]<-p[2]*(f.pos[t] *(1-psi)*(1-phi[t])  +  f.pos[t] * psi*(1-v)*(1-phi[t]))  #changes


#row 2
A[2,1]<-p[1]*(1-m)*(1-phi[t])
#row 3
A[3,2]<-p[2]*(1-phi[t])
A[3,3]<-p[2]*(1-phi[t])

#row4 
A[4,3]<-p[2]*f.neg[t]*phi[t]
A[4,6]<-p[2]*f.con[t]*(v+phi[t]-v*phi[t])
A[4,7]<-p[2]*(f.pos[t]*psi*(v+phi[t]-v*phi[t]) + f.pos[t]*(1-psi)*phi[t])


#row5
A[5,1]<-p[1]*phi[t]*(1-m)
A[5,4]<-p[1]*(1-m)

#row6
A[6,2]<-p[2]*phi[t]
A[6,3]<-p[2]*phi[t]
A[6,5]<-p[2]
A[6,7] <- 0

#row 7
A[7,6]<-p[2]
A[7,7]<-p[2]


#row 8
A[8,1]<-p[1]*m
A[8,4]<-p[1]*m
A[8,8]<-p[3]
#Do eigen analysis for stable age distribution
	ea=eigen.analysis(A)
 	
 	
	stage[j,]=ea$stable.stage
	adult.fem[j] = sum(stage[j,c(3,6,7)])
	yearling.fem[j] = sum(stage[j,c(2,5)])
	adult.male[j] = stage[j,8]
	calf[j] = sum(stage[j,c(1,4)])

lambda.diff[j] = lambdaH[j] - lambdaI[j]
	
}





lambda_healthy=list(lambda=lambdaH,adult.fem=adult.fem,yearling.fem=yearling.fem,calf=calf,adult.male=adult.male)

lambda_infected=list(lambda=lambdaI,lambda.diff=lambda.diff,adult.fem=adult.femI,yearling.fem=yearling.femI,calf=calfI,adult.male=adult.maleI)

save(lambda_healthy, file ="lambda_healthy.Rdata")

save(lambda_infected, file = "lambda_infected.Rdata")
#
#